We present a methodological update to the Pythonic Black-box Electronic Structure Tool (PyBEST), a fully-fledged modern electronic structure software package. Among other things, these include the extensions to pair Coupled Cluster Doubles (pCCD) open-shell electronic structures via the ionization potential and double ionization potential equation of motion formalism (IP-EOM-pCCD and DIP-EOM-pCCD), tailored Coupled Cluster methods, EOM extensions on top of the linearized Coupled Cluster correction on top of pCCD, the configuration interaction correction on top of pCCD, and the static pCCD-in-DFT embedding. Furthermore, we have increased the performance of PyBEST by optimizing its underlying tensor contraction engine and providing GPU support via the CuPy library. New version program summaryProgram Title: PyBEST: Pythonic Black-box Electronic Structure ToolCPC Library link to program files:https://doi.org/10.17632/xf9kb7yfwr.2Developer's repository link:https://pypi.org/project/pybest/ and https://zenodo.org/records/10069179Licensing provisions: GPLv3Programming language: Python3.9–Python3.12, C++Journal reference of previous version: Computer Physics Communication 264 (2021) 107933.Does the new version supersede the previous version?: YesReasons for the new version: The new version includes significant speedup and memory optimizations. We provide an interface to CUDA GPU via the Python CuPy library. New methodological features are the ionization potential/double ionization potential equation of motion pair coupled cluster doubles (IP/DIP-EOM-pCCD), configuration interaction on top of pCCD (pCCD-CI), tailored CC methods including DMRG-tCC, density matrices from LCC(S)D and pCCD-LCC(S)D, a static WFT-in-DFT embedding and an external charge module. We also introduced support for frozen core orbitals in correlated methods, dumping orbitals in the cube format, printing of Hartree–Fock (HF) orbital energies, the one-dimensional contact potential model Hamiltonian, and wrappers to calculate properties such as dipole and quadrupole moments. To do that, we redesigned and optimized the class structure of several PyBEST's core sub-packages, such as the linear algebra, geminals, and occupation model classes.Summary of revisions:1.GPU support (via the CuPy library [1]) for selected contractions in PyBEST's [2] tensor contraction engine [3]. It allows us to obtain a significant speedup (up to a factor of 15) for bottleneck contractions in standard coupled cluster and post-pCCD [4] implementations in PyBEST using the Cholesky linear algebra factory. Overall, a speedup of approximately a factor 4 is achieved [3]. We provide a batch-wise algorithm, which alleviates the limitations of a single Video Memory (VRAM) and enables large-scale quantum chemical calculation even on one GPU card [3].2.Several new quantum chemistry methods have been implemented in PyBEST. These include:•EOM-LCCD and EOM-LCCSD methods with an RHF/pCCD [5–9] reference function for the calculation of electronic excitation energies [10–12]. The implementation allows for targeting singlet-singlet excitation energies, including doubly excited states.•IP-EOM-pCCD and DIP-EOM-pCCD methods [13] that allow us to compute open-shell electronic structures on top of a pCCD reference function. Our implementation supports singlet, doublet, triplet, quartet, and quintet spin states. The user can choose between a low-spin and different high-spin versions of these methods.•A CI-based dynamic energy correction on top of RHF and pCCD: pCCD-CID and pCCD-CISD only support a spin-free implementation, while for CIS, CID, and CISD, the user can choose between the Slater Determinant (SD) and Configuration State Function (CSF) implementations. We feature various Davidson-type corrections on top of a CI wave function [14].•tailored CC-based methods [15] that use a pCCD wave function as the fixed reference function (so-called frozen-pair coupled cluster methods [16,15]) or the external amplitudes are provided from different wave function ansätze, e.g., from an MPS representation or any external file.•Density matrices from LCC(S)D and pCCD-LCC(S)D methods [17]. These are used to determine orbital entanglement and correlation measures [18–21, 17] and other method-dependent properties.•Wrappers for calculating nuclear and total dipole and quadrupole moments using integrals provided in the Libint library.•A static embedding module, where all wave-function-based methods provided in PyBEST can be embedded in the potential generated from external software [22].•Introduction of a universal and more general way to define core orbitals in post-HF calculations.3.New Hamiltonians:•Support for 1-D contact interactions (Hˆii′=Tˆii′+Vˆii′) in an arbitrary external potential (Vˆii′=δii′V(xi)) [23]. Tˆii′ represents the 1-D kinetic energy, and xi=iΔx is a grid point, where i can take values of 0,±1,±2,…. Thus, this choice of the Hamiltonian leads to a uniform grid in the coordinates. It corresponds to using a Fourier basis set.•Augmentation of the quantum chemical Hamiltonian with external point charges [22]. The following one-electron integrals have been added:(a)integrals for the interaction of point charges and electrons(b)integrals for the interaction between point charges and nucleiHˆ=Hˆel+∑i∫ρel(r)qidr+∑i∫ρnuc(r)qidr, where Hˆel is the electronic Hamiltonian of the molecular system, qi is the i-th external point charge, ρel and ρnuc denote electronic and nuclear densities of the molecular structure under consideration.(c)the point charge counterpart of pVp integrals. For the relativistic Douglas–Kroll–Hess Hamiltonian in PyBEST, the picture-change error correction is also provided [24].4.Code optimization:•Replacement of slicing and C++-based routines with numpy.tensordot, numpy.einsum, and opt_einsum, which leads, for instance, to a speedup of about a factor of 50 when constructing the Fock matrix.•Reduced memory requirements in post-pCCD calculations. We allow for dumping effective Hamiltonian objects in all CC sub-packages to disk if the intermediates are larger than o2v2 to reduce RAM consumption in large-scale calculations [25], where o is the number of occupied orbitals and v the number of virtual orbitals with respect to some reference determinant.•New flavors of the 4-index transformation based on opt_einsum and a naive einsum implementation [26] and replacement of former slicing operations.•Memory-optimized Cholesky routines throughout the codebase, which reduce the memory peak even by a factor of 2.•New select features for PyBEST's contract function in its tensor contraction engine. It defaults to CuPy for specific bottleneck operations.5.Code restructuring:•A more optimal structure for the linear algebra sub-package that avoids circular imports.•All .h5 format interactions are now done purely with h5py.•Support for dumping the occupation model instance to a .h5 file.•Rewrite of occupation model class to support future updates and extensions (automatic definition of frozen core orbitals and the possibility to introduce an active orbital space).•A new structure of the geminal module to allow for future support of other geminal models.•Change of the contraction order for tensor contractions (cupy → tensordot → einsum/opt_einsum: [Display omitted] Nature of problem: The idea behind the PyBEST software package is the efficient and reliable modeling of electronic structures and properties of small and large-scale quantum systems at the interface between chemistry and physics using geminal-based, and specifically pCCD-based, methods. Such models represent an efficient and compact representation of the wave function. The aim of PyBEST is focused on reliable predictions of organic electronic materials [27,25] and actinide-containing complexes and their properties [28–30].Solution method: The modular design and general tensor contraction engine implementation in PyBEST allows for a convenient and prompt code acceleration by drop-in-replacements of modern Python-based packages. Combined with a series of unconventional (and conventional) electronic structure models based on the pCCD ansatz to solve the electronic Schrödinger equation, PyBEST represents a modern alternative to the existing quantum chemistry and physics software packages. With PyBEST, not only can we model ground- and exited-state properties, but also provide a better understanding of orbital interaction through orbital-based quantum-information measures.Additional comments including restrictions and unusual features:PyBEST features unique electronic structure methods based on orbital-optimized pCCD and their interpretation in terms of orbital entanglement and correlation that are not available in any other software package. While pCCD-based methods are limited to closed-shell systems, open-shell configurations can be accessed via the IP/DIP-EOM extensions.